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Abstract 

Using harmonic and anharmonic force constants extracted from density-functional calculations 
within a supercell, we have developed a relatively simple but general method to compute thermo- 
dynamic and thermal properties of any crystal. First, from the harmonic, cubic, and quartic force 
constants we construct a force field for molecular dynamics (MD). It is exact in the limit of small 
atomic displacements and thus does not suffer from inaccuracies inherent in semi-empirical poten- 
tials such as Stillinger- Weber's. By using the Green-Kubo (GK) formula and molecular dynamics 
simulations, we extract the bulk thermal conductivity. This method is accurate at high tempera- 
tures where three-phonon processes need to be included to higher orders, but may suffer from size 
scaling issues. Next, we use perturbation theory (Fermi Golden rule) to extract the phonon lifetimes 
and compute the thermal conductivity k from the relaxation time approximation. This method is 
valid at most temperatures, but will overestimate k at very high temperatures, where higher order 
processes neglected in our calculations, also contribute. As a test, these methods are applied to 
bulk crystalline silicon, and the results are compared and differences discussed in more detail. The 
presented methodology paves the way for a systematic approach to model heat transport in solids 
using multiscale modeling, in which the relaxation time due to anharmonic 3-phonon processes is 
calculated quantitatively, in addition to the usual harmonic properties such as phonon frequen- 
cies and group velocities. It also allows the construction of accurate bulk interatomic potentials 
database. 

PACS numbers: 63.20.-e,63.20.dk,63.20.kg,61.50.Ah 



I. INTRODUCTION 

Classical molecular dynamics (MD) simulations use either semi-empirical potentials such 
as Stillinger- Weber (SW)[l|, Abell-Tersoff-Brenner|2| or other type of force fields where 
the potential energy is an analytical function of the atomic positions, or first-principles 
potentials calculated typically using density-functional methods based on either the Born- 
Oppenheimer[31] or the Car-Parrinelloj^ dynamics. The former are fast to compute but suffer 
from inaccuracies, while the latter are accurate but time-consuming to compute. Due to 
recent interest in thermal transport in semiconductor materials having good thermoelectric 
properties and the topic of microelectronics thermal management in general, there have 
been many calculations of the lattice thermal conductivity of materials using the Green- 
Kubo (GK) formula 5|, IS]. This formula relates the thermal conductivity, through the use of 
the fluctuation-dissipation theorem, to the time-integral of the heat current autocorrelation 
function. The latter is calculated from an MD simulation, and the ensemble average is 
usually replaced by a time average. Semi empirical potentials such as SW are usually used 
to perform the MD simulation for a system such as Si. As the thermal conductivity of 
a perfect crystal is mainly due to anharmonic three-phonon processes, directly related to 
the third-derivatives of the potential energy with respect to atomic displacements, and the 
latter is generally not fitted or considered in the design of the semi-empirical potentials, 
there is really no good reason to expect an accurate value for the thermal conductivity 
calculated from a GK-MD simulation. In the case of Si, when using the SW potential, 



however, for some reason [7|, relatively good agreement is found between the experiment and 
the simulation results, even for a relatively small sup er cell |8Nll|. The latter fact is also 
cause for concern, because, as we will show in the following, a small supercell limits the 
number of long-wavelength phonons which carry a large portion of the heat in a material. 
The lucky agreement can be attributed to a cancellation due to two different effects, which 



will be discussed in sectioi JITTBl Sellan et al. [Uj have presented a discussion on convergence 



issues, mostly with respect to non-equilibrium MD simulations, and we will also discuss the 
scaling issue with respect to GK-MD and lattice dynamics (LD) in this paper. 

More accurate calculations of the thermal conductivity, based on the full solution of 
Boltzmann transport equation have shown that the thermal conductivity of Si using the 
SW potential is about 4-times larger than the experiment, while using the Tersoff or EIDP 



potentials produce results that are about twice larger than the experimental values 12 1. 



In a similar work, using the environment-dependent interaction potential (EDIP), Pascual- 



Gutierrez et al. [l3l| also find thermal conductivity of bulk Si from MD in good agreement with 
experiments. In a subsequent work[14], the same group computed the thermal conductivity 
using the lattice dynamics (LD) theory based on the same EDIP potential, similar to the 



work of Broido et al. [12|. They obtain good agreement with experiments whereas Broido 
et al. do not. The reason for this discrepancy, as they also mention in their paper, is 
unclear. Opinions on the accuracy of semiempirical potentials seem to differ as there have 



been other reports 15| where SW is found to overestimate the thermal conductivity by 70%. 
As we will show in this paper, some of these potentials might not be completely reliable 
for the calculation of the thermal transport properties for the simple reason that they were 
not fitted or constructed to have the correct third derivatives, which are responsible for 
the thermal resistivity of a material. In fact even their harmonic force constants produce 
phonon dispersions and elastic constants which differ from experiments by 10 up to 40%. 
Furthermore such potentials exist and have been thouroughly tested for only a very small 
number of pure crystalline solids. 

In a tour de force work, Broido et al.., later, used the density- functional perturbation the- 
ory (DFPT) formalism in order to calculate the phonon scattering rates from first-principles 
DFT calculations and were able to successfully reproduce the thermal conductivity of bulk 
Si and Ge 16|. Their approach, which was very accurate, included the calculation of all the 
cubic force constants up to 4 lattice parameters away, and the complete iterative solution of 
the Boltzmann transport equation. 

We recently developed a methodology to extract second, third and fourth derivatives 



of the potential energy from first-principles calculations ITJ], and showed that the phonon 
dispersion relation in Si can be well-reproduced. In this paper, we pursue this work further 
and use these derivatives to construct a force field in order to explore the results from MD 
simulations and perturbation calculation to calculate the thermal conductivity of bulk Si. 
We should mention that our approach, even though very similar in essence to that of Broido 
et al., is simpler in the sense that it limits the range of the force constants (FCs) to a few 
neighbors (5 for harmonic and 1 for cubic in the present case study of Si, in contrast to 



more than 20 for harmonic and 10 for cubic in the work of Broido et al. 16|]). For the sake of 



physical correctness, however, we enforce the translational, rotational and Huang invariances 



on the extracted force constants. So the latter are not exactly equal to the ones obtained 
from DFPT or any finite difference calculation of the forces, but make the calculation load 
much lighter than if one had to include so many neighbors. On the other hand a DFPT 
calculation, if restricted to few neighbors, should enforce all these invariances. Usually, only 
the translational ones, also known as the "Acoustic Sum Rule" (ASR) are enforced in some 



of the standard DFT codes such as Quantum Espresso |18[|. 

In what follows we briefly review the methodology to extract force constants from first- 
principles density-functional theory calculations (FP-DFT). The formalism for the molecular 
dynamics and Green-Kubo calculations of n are explained in section fillip . This will be 
followed by the lattice dynamics approach detailed in section (lIVp . Results for Si will be 



shown and discussed in section (jV|), followed by conclusions. 

II. EXTRACTION OF FORCE CONSTANTS FROM FP-DFT CALCULATIONS 

We construct the potential energy V for the MD simulation as a Taylor expansion up to 
fourth order in the atomic displacement Ui of atom i about its equilibrium position: 



V = Vo + J2 ^i^i + ^ II ^ij ^^^0 (1) 

i ■ ij 

+ T]l^ ^ijk ^i^i^k + jy Z^ Xijkl UiUjUkUl 
ijk ' ijkl 

= E(e^ - l^^^f) (2) 

The last equation defines the on-site energy Cj . If the displacements Ui are around the 
equilibrium position, and this is usually the case, IIj = 0. The coefficients $, \E' and x i^i 
the expansion are called the harmonic, cubic, and quartic force constants respectively, and 
satisfy certain symmetry constraints. Namely they must be invariant under interchange of 
the indices, uniform translations and rotations of the atoms, in addition to invariance under 



symmetry operations of the crystal. The detai 



imposed can be found in our previous work 17|. 



s of the needed constraints and how they are 



To get these numbers, we consider one or several supercells in which atoms are in their 
equilibrium position. One, two or three neighboring atoms are moved simultaneously by a 
small amount, typically about O.OlA along the cartesian directions. Consideration of crystal 



symmetry usually reduces the needed displacements. For instance in a cubic-based crystal, 
like silicon, where the two atoms in the primitive cell are equivalent, one only needs to 
move one Si atom along the x direction. This is sufficient to extract all harmonic force 
constants if the supercell size is large enough. The latter size is chosen depending on the 
available computational power and the considered range of force constants. To get three- 
and four-body interaction terms, one needs to move two and three atoms at a time and 
record the forces on all atoms in the supercell. It is advantageous to record forces for 
atomic displacements in two opposite directions as there would be cancellation of the cubic 
contributions and this would make the calculation of the harmonic FCs much more accurate 
(up to order u"^). Thus one would obtain a large set of force- displacement relations computed 
from a FP-DFT code. Together with the invariance constraints, an overcomplete linear set of 
equations on all the force constants will be formed. A singular- value decomposition algorithm 
is then used to solve this linear overcomplete set. We find that usually the violation of the 
invariance relations is of the order of 10~^ times the FC itself. This however requires very 
accurate evaluation of the forces, meaning that they should have converged with respect 
to the cutoff energy and number of k-p oints to within at least 4 significant figures, if not 
more! Our experience on graphene 19| and silicon (present work and reference [ll7|] has 
shown that the harmonic force constants are usually reproduced quite accurately. Higher 
order FCs have less accuracy as their contribution shows up in the second or third or fourth 
significant figures of the forces. The main approximation is in cutting off the range of the 
interactions, which will lead to inaccuracies in some of the Gruneisen parameters as we will 
shortly see. 

III. MOLECULAR DYNAMICS AND THE GREEN-KUBO FORMALISM 

Typically a supercell is constructed with periodic boundary conditions, and an MD sim- 
ulation is performed over a long enough time steps in order to reach thermal equilibrium, 
followed by a long (N,V,E) simulation in order to collect data on J for later statistical 
processing, i.e time and ensemble averaging of its autocorrelation. 

Based on the potential displayed in Eq. ([2]), one can extract the expression for the force, 
required in the MD simulations: 



The heat current is defined in the discrete (atomic) case of a lattice, where there is no 
convection, as: 

^" = E ^r = E ^^ = E(^r - ^")(^. ■ ^) (4) 

where Vi is the velocity of particle i and e,, as defined in Eq. ([2]) is the local energy of atom 
i. Using our expansion, it can be expressed as a function of the force constants as follows 
(as we expand around the equilibrium position, we assume IIj = 0) : 



+ Y^Y^XijkiUjUkUi] (5) 

ijkl 

This definition leads to the following form of the local heat current: 



1 

+ jY.XijkiUkUi] (6) 

^ kl 

^''=Vk^L <J°W'V)>dt (7) 

where a, /3 = x, y, z, and < A> denotes the equilibrium average of the observable A, which 
in the classical case can be replaced by its time average provided the time is long enough to 
satisfy ergodicity. The ensemble averaging is necessary as long as different runs starting 
with different initial conditions lead to different integrated autocorrelation functions. The 
true current autocorrelations decay quite fast. In a single MD run, however, this decay is 
not observed. Instead one observes a decay in the amplitude followed by random oscillations 
about zero. As finite systems are usually not ergodic, an ensemble average, over the random 



initial conditions is also needed to correctly simulate the equilibrium average required in 
the Green-Kubo formula. In this case, the ensemble-averaged autocorrelation function will 
decay smoothly to zero with time. One will then see that the decay is indeed relatively short, 
because the long-time tails get cancelled after ensemble averaging. The advantage of the 
ensemble averaging, in addition of course to the usual time averaging, is that one samples 
the phase space more randomly, and generates uncorrelated sets of pairs J°'{0)J^{t) for a 
given time difference t. In this case, the mean has the convergence properties of gaussian- 
distributed variables and the error decays to zero as the inverse square root of the number 
of initial conditions. 

In a numerical simulation, the GK formula should be replaced by: 

1 1 ^= At 

[ E^E JtipU^^P + t)] (8) 

t=0 p=l 

where T is the total simulation time of each run. At is the time difference between two 
successive data points, t and p are integers labeling time, and A^ens is the number of generated 
initial conditions, each labeled by i. One important comment is in order here, and that is 
the use of 1/T in the denominator instead of the more intuitive 1/(T — t + 1) which is the 
actual number of terms in the last sum. One can show that the choice in Eq. (jHj) provides 
an unbiased estimator of the autocorrelation [20|]. Some previous works in the GK method 
such as reference 9J] have used the "biased" formula, which would be fine as long as the 
total simulation time is much larger than the largest time t used for the integration. One 
simple way to become convinced that Eq. ([8]) is the correct way, is to consider large times 
t near the maximum simulation time T. As J{p) and J{p + 1) are both fluctuating random 
numbers, and there are not enough terms in the sum to make the average go to zero, the 
time average will become large again and tend to aj =< J(0) J(0) > instead of decaying to 
zero as t approaches T. A division by the small number T — t would not solve the problem, 
whereas the division by T would make this very small as it should. 

The simulation time T being usually quite long, the statistical error in the time averging 
is usually small, and we estimate the error bars in our data from the ensemble averaging: 
If C{t) is the ensemble-averaged autocorrelation function, its error bar AC is evaluated as: 

8 



AC{t) = [Etrmt) - c(t))2]V7iVe„, 

The magnitude of this error bar and the required accuracy in the results determine how 
many ensembles are needed for a proper estimation of thermal conductivity. 

A. How many MD steps are necessary? 

For a given supercell size, there is a discrete number of phonon modes which can propa- 
gate and get scattered in the system. The largest wavelength consistent with the periodic 
boundary conditions would be the supercell length. To this, one can correspond a small- 
est phonon wavenumber or frequency allowed in the simulation: Ucut = ^nc/L. The total 
simulation time should be large enough so that all phonon modes can get scattered a few 
times within the simulation period. Largest relaxation times belong to acoustic modes, and 
usually decay as the inverse square of the phonon frequency. Knowing the smallest allowed 
frequency Ucut due to the finite size of the system, one can estimate the corresponding 
phonon lifetime (using Klemens' formula displayed in Eq. (Q for instance). The total MD 
simulation time should therefore be a few times larger than the largest phonon lifetime so 
that scattering events of long wavelength phonon modes can properly be sampled during 
the MD run. As an example, we can consider Si system in a cubic 10x10x10 supercell of 
8000 atoms. In this case L = 10 x 5AA = 5.4 nm corresponding to kc = 27r/10a which is 
a fifth of the F — )■ X line. The lowest frequency mode is therefore about Ucut = 1 THz. 
As can be seen in Fig. ([5]) below, the normal and umklapp lifetimes at this frequency are 
3000 ps and 10000 ps respectively, leading to a total lifetime of 2300 ps. So in order to sam- 
ple such rare events, one needs to calculate the autocorrelation over at least 20 ns, which 
means the MD simulation should be run for at least the same amount of time if not longer. 
This could be computationally prohibitive. If the runs are made with fewer MD steps long 
wavelength phonons will not be relaxed and the autocorrelation would not tend to zero. 
Another consequence of this remark is that for large supercells, as relaxation times scale as 
l/cj^^j oc L^, longer simulation times proportional to the square of the the supercell size 
would be required. 



B. Size scaling 

As mentioned, the choice of a finite supercell conies with the cost of discretizing the 
phonon modes and supressing the phonons of wavelengths longer than the supercell length. 
The neglected contribution maybe estimated as follows: the anharmonic lifetime of acoustic 
modes maybe approximated by the Klemens' formula |30| 



1 _ 2 SfcfiT u;2 

^Klemens 'kX t\t 2 . .max 



2 --»- -fcA /qA 

ikX /\/f „,2 , .max \'^J 

' kX ^^^ ^kX ^x 

where u^""^ is the largest frequency of the branch A, 'jkx the mode Gruneisen parameter, 
Ukx the frequency, and Vkx the group velocity associated with the mode kX. Therefore long 
wavelength phonons will have a large relaxation time and can considerably contribute to 
the thermal conductivity. Assuming this form in the relaxation time approximation to the 
thermal conductivity, and using Eq. (ITTll . we can write the thermal conductivity as a sum 
over contributions of phonons of different frequencies: 

/•oo 1 

n= / -T{ujy{Lo)C^iLo)DOS{Lo)duj 
Jo 3 

In 3D, since the density of states (DOS) is quadratic in frequency, the contribution of 
long wavelength acoustic phonons would be linear in the cutoff frequency Ucut = 27ic/L: 

k{L) = /t(oo) - r™' -t{uj)v\u)C^{u)DOS{uj) du 
Jo 3 

For low frequencies Cv{u) = kB[f3hco / 2smh.{f3hijj / 2)]"^ ~ ks and v{u) ~ c so that 

r^cut \ 

k{L) = k{oo) —A -—DOS{uj) du = k{oo) — DoOcut 

Jo LO^ 

= k(cx)) -E- = /t(cx)) - F/yt (10) 

where A,D,E and F are constants which do not depend on the size, A^ is the mean free path 
(MFP) associated with the cutoff frequency Ucut = 2ttc/L. This gives us a way to deduce 
how the thermal conductivity of a finite size sample scales with the supercell length or the 
cutoff frequency Ucut or MFP, Ac. We should note that a different scaling law {1/k{L] 



1/k{oo) + C/L) was also proposed and used by Sellan et a/.[ll[, and also Turney et al.[21\ 
when they want to extrapolate their thermal conductivity data to infinite size. The argument 
they used to deduce it however was based on the Matthiesen's rule, stating that the bulk 



10 



resistivity is obtained by adding the finite size resistivity to the one obtained from L/v taken 
as relaxation time. 

There is another additional problem with finite size MD simulations. Even though mo- 
mentum is still conserved in a 3-phonon process, because the modes are discrete in a 
finite supercell, energy conservation will not always be possible, unless the energy difference 
uj — uji — UJ2 <T where F is on the order of the sum of inverse lifetimes of the three con- 
sidered phonons. If this relation is not satisfied, the considered 3-phonon scattering will not 
take place in a finite supercell, and this will lead to an overestimation of the lifetime of the 
phonons, and thus, of the thermal conductivity. 

These competing effects, namely an overestimation of k due to limited phase space for 
energy conservation and an underestimation due to cutoff of low frequency acoustic modes, 
may lead to a magical cancellation, resulting in thermal conductivities in good agreement 
with experiments even for moderate supercell size. This error cancellation will likely affect 
the temperature-dependence of k: at higher temperatures the discreteness error is reduced 
as r increases linearly with T. The frequency cutoff error, however, will not be affected 
by high temperatures. Consequently, as T is increased the thermal conductivity of a finite 
sample will decrease faster than 1/T with temperature. This has been observed in the work 
of Volz and Chen[8:]. It can also be verified by introducing other scattering events such as 
isotope or defect scattering leading to larger V values. In such cases, the discreteness of 
modes will have little effect, and the simulated k will be less than the exact one, due to 
the cutoff of long MFP phonons effect. As a result, in a system where due to disorder or 
high temperatures scattering rates are high, GK-MD simulations will typically require larger 
supercells to converge. 

The correct way of estimating i^{T) is to do a proper size scaling at each temperature by 
plotting k(T, L) versus 1/L and linearly extrapolating to 1/L — )• 0. 

IV. THE LATTICE DYNAMICS APPROACH 

Using the extracted for constants, one can form the dynamical matrix of the crystal using 
its primitive cell data: 
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where R is a translation vector of the crystal, r refers to an atom in the primitive cell, and 
a, P are cartesian components x, y, z. Such sums of the force constants over the translation 
vectors of the primitive lattice are usually short-ranged and fast to compute, except if 
Coulomb interactions are involved, in which case the sum is evaluated using the Ewald 
method. 

Diagonalizing this matrix, one can find the phonon spectrum and the normal modes as 
its eigenvectors: 

ED:^'if^)<^ik)=^LeT{k) (12) 

r'/3 

where A labels a phonon band (or branch), and k refers to a point in the first Brillouin zone 
(FBZ). Using these eigenvectors and eigenvalues, and from perturbation theory, one can 



calculate the phonon 



self-energy defined as[2J- 



ineshifts and lifetimes as the real and imaginary parts of the 3-phonon 



m- 



S(gA,u;) = --^ ^ \V{qX,l,2)\' x 
^1 + ni -t- na) (n2 - rii) 



(13) 



ui + U2 + euc ui — U2 + euc 

where Uc = oo — irj , (r/ ^ O"*") is a small infinitesimal number, which in practice is taken to 
be finite for a given k-mesh size, n is the equilibrium Bose-Einstein distribution function, 
and 1 and 2 refer to modes (qiXi) and (g2A2)- The 3-phonon matrix element V, expressed 
as a function of the cubic force constants \E', is given by: 

e^^^"^^+'^"^^^er(9)eir(9i)eir(g2) . . 

The calculation of the self-energy would require a double sum over the q-points (labeled 
above by 1 and 2) in the FBZ . Due to the conservation of momentum, however, only 
terms with q + qi + q2 = G, with G being a reciprocal lattice vector, should be included 
in the above sum. In practice, therefore, this involves only a single summation. To get 
the phonon dispersion and lifetimes due to 3-phonon scattering terms, one needs to solve 
E = ujkx + Tj'{k\,E) where S' is the real part of the self-energy (S = S' + iS") . This 
equation needs to be solved iteratively. Since the shift is usually small, to leading order, 
one can use E = Ukx + T,'{kX,Ukx) i-e. one uses the on-shell frequency as argument of the 

12 



self-energy. The same approximation will be used for the imaginary part giving the inverse 
lifetimes. The corresponding phonon lifetime will be given by Tk\ = 1/21]" {kX,uJkx)- In 
the evaluation of the imaginary part S" , one encounters Dirac delta functions reflecting the 
conservation of energy in the three-phonon process: u^x = ooi ±^2. In effect, from Eq. flT^ . 
it can be noticed that the delta function is substituted by a Lorentzian function of width 77. 
The latter depends on the choice of the k-point mesh in the FBZ. A small value for 77 can be 
used for a fine mesh, while a coarse mesh requires larger values of rj. Typically rj is chosen 
to be of the order of energy spacing in the joint density of states (JDOS) so that the latter 
is a smooth function of the frequency and does not display any oscillations with sharp peaks 
which would appear if the width is too small. 

JDOS{uj) = — — ^ 6{uj — uji — U2) + S{u — uji + U2) (15) 

The anharmonicity can be characterized by the Gruneisen parameters (GP). The force 
constant GP is defined as 'j^ = -d In 0/2 d In \^ where V is the volume. The mode GP 
is defined as: •ykx = -d In cokx/ d In l^ where Ukx is the phonon frequency evaluated at the 
point k and band index A. It gives the relative decrease in the phonon frequency as the 
volume is increased by 1%. From the Taylor expansion of the harmonic force constants in 
terms of the volume or the lattice parameter, one can calculate such change. 



7.A = -7-^E^o"} 



«2 
RlTi,R2T2 



e 



ik-{R2-Ri) 



X X^^er^{-k)er'{k) (16) 

where Xji^- is the equilibrium atomic position of atom type r in the primitive cell labeled by 
the translation vector R. 

Finally, the thermal conductivity is calculated within the relaxation time approximation 
(RTA), which leads to the following well-known expression for the thermal conductivity: 

'^ = ^?T^ 51 ^L^fcA ^ujkx dukx/dT (17) 

where ft is the volume of the unit cell. The relaxation time t^x in this expression represents 
the time after which a phonon in mode kX reaches equilibrium on the average, and depends on 
the scattering processes involved. In a pure bulk sample, the only source of phonon scattering 
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is anharinonicity dominated usually by tliree-phonon processes. Using perturbation theory 
or the well-known Fermi Golden rule (FGR), one can derive the expression of the relaxation 



24 



27] . It can be shown that to a good 



time as a function of the cubic force constants 
approximation, it is given by 

^'^ "^ 23[S(gA,M] ^^^^ 

In what follows, we have disregarded the boundary scattering term, which is responsible 
for the low-temperature behavior of k. In such case, k is expected to saturate to a finite 
value at low enough temperatures. The reason for this saturation can be understood if one 
assumes the low-frequency limit of the DOS and the relaxation times similar to Eq. [TUl 
Considering that in w — )■ limit we have: DOS\{uj) — )• u"^ /2t:'^c\ and C^{(jj) = ^^(x/shx)^ 
(with X = (3hu/2), and the relaxation time can be written as t{uj) — )■ hu^/u'^kBT, the 
integral defining the thermal conductivity can be transformed, in the low temperature limit, 
to: 



tt^ca Jo shx Qcx 

The constant Uo that appears in the low energy limit of the relaxation time as well as the 
speeds of sound c\ determine the saturated value of the thermal conductivity. So when 
only 3-phonon scattering processes are included, the thermal conductivity would tend to 
ksOjl/Gcx as T goes to 0. 

Finally, in our numerical calculations where the integral in the FBZ has been approx- 
imated by a sum over a discrete set of k-points, the low-frequency region is not properly 
sampled and we observe a decay to zero at low T, and therefore have not reported the 
unreliable low-temperature data in this work. 

V. RESULTS AND DISCUSSIONS 

A. Validation of force constants 

First-principles calculations were done using the PWSCF code of the Quantum Espresso 
package [18!] A set of force-displacement data were calculated using 2x2x2 supercell of 64 
Si atoms. The set of force-displacements data, along with the symmetry constraints, form 
an over-complete linear set of equations needed to determine the potential derivatives. We 



use the local density approximation (LDA) of Perdew and Zunger 22| with a cutoff energy 
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TABLE I: Typical deviations in the SW and Taylor expansion (present model) force fields compared 
to true FP-DFT forces. They are obtained by moving all 64 atoms in the supercell in a random 
direction by 0.1 and 0.2 A respectively. 



Amplitude(A) 


ct(SW) 


cr(Present) 


0.1 
0.2 


0.35 

0.28 


0.05 
0.08 



of 40 Ryd and 10 k-points in the irreducible Brillouin zone of the cubic supercell. The 
range of different ranks of force constants can be chosen by the user. We have set the 
range of harmonic forces constants (FCs) to 5 nearest neighbor shells, and that of the cubic 
and quartic force constants to the first neighbor shell only. This results in 17, 5 and 14 
independent harmonic, cubic and quartic FCs respectively. The corresponding number of 
terms in the Taylor expansion of the potential energy are, however, equal to 1500, 1146 and 
7980 respectively. This is why the ranges were restricted to 5, 1 and 1 nearest neighbor shells 
in order to limit the computational time to a reasonable amount. Note that despite the large 
number of terms to be computed, arithmetic operations are only limited to additions and 
multiplications. 

In Fig. ([T]), we show the change in the total energy as an atom in the supercell is moved 
along the [100], [110] and [111] directions respectively. Resutls from DFT calculations are 
compared against our developed force field including the harmonic, harmonic+cubic, and 
harmonic+cubic+quartic terms of the Taylor expansion. For the sake of comparison, we have 
also plotted the same energy change ontained from the Stillinger- Weber (SW) potential|l[, 
which is widely used in MD simulations of Si systems. 

To further assess the accuracy of the force field, we have also moved all the atoms in 
the supercell in different random directions by a small amount of magnitudes 0.1 and 0.2 
A respectively, and compared the average force of our model and the SW potential to the 
FP-DFT one. The deviation is charaterized by: 



E/pmodel _ pDFT\ 

a[uioae\.) = — 



y^ pDFT2 



(19) 



The results for the parameter a are summarized in table ([T]). 

We can notice that this type of error estimate would also include contributions from 
many-body forces, and is a more stringent test on the force field. The errors from the 
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FIG. 1: (Color online) Total energy as an atom is moved in the [100] (left), [110] (middle) and 
[111] (right) directions. DFT results are compared with the force field and the Stillinger- Weber 
potential. MD234 refers to the force field in which all harmonic, cubic and quartic terms are 
included, while MD2 refers only to the harmonic force field, etc. 

present model are consistently about 4 to 5 times smaller that the SW potential. 

In the following we follow two paths to compute the thermal conductivity. The first is to 
use the Green-Kubo formula, by using the results from an MD simulation: 

B. Thermal conductivity from MD-GK 

As previously mentioned, there will be large fluctuations in the current autocorrelation 
function versus time from one run to the next, and therefore an averaging over several initial 
conditions is necessary to produce a reliable plot. In Fig. ([2]), we have plotted such ensemble 
average for a 10x10x10 supercell containing 8000 atoms. The error bars are mainly due to 
the ensemble averaging, and those related to the time averaging are small as the number of 
MD time steps are quite large. 

We can also see in this figure the cumulative integral of the ensemble-averaged autocor- 
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relation function. The same calculation was performed in a 7x7x7 supercell of 2744 atoms, 
where the averaging was over 99 runs with different initial conditions. Due to its larger size, 
there are smaller fluctuations in the average current per atom in the 10x10x10 supercell, and 
we only used 27 initial conditions for this supercell. Since from each MD run one can really 
extract three autocorrelation functions k^x, Kyy and Kzz, which are equal by cubic symmetry, 
we also averaged over the 3 directions. In this sense, the above mentioned numbers should 
be multiplied by 3. 



Autocorrelation and its cumulative integral (10x10x1 

80 

Integrated autocorrelation 
Autocorrelation 




300 400 
MD time {ps) 



FIG. 2: (Color online) Plot of the ensemble-averaged (over 27 initial conditions) heat current 
autocorrelation as a function of time, and its integral for the 10X10X10 supercell. Vertical units 
for the integrated autocorrelation are in W/mK, and the autocorrelation (blue dots) has been 
multiplied by a constant to be on scale. 



The error bars are determined by the large fluctuations in the integrated autocorrelations 
divided by the square root of the number of ensembles. The error bar due to the time average 
is usually much smaller if MD simulations are run for a long enough time. 

The results for two different supercell sizes are summarized in table (ITll) as compared 
with the experimental data of Slack et a/. [28]. One can notice an underestimation of the 
experimental data, which is reduced as the supercell size is increased. To get the correct 
value in the thermodynamic limit, one needs to extrapolate these results to infinite size. 

There are a few competing effects which can explain this discrepancy: the most important 
one is size effect, which as was just explained, underestimates k. Similarly, the larger value 
of the Gruneisen parameter for the acoustic modes in our model will produce a smaller 
relaxation time (see the Klemens formula in Eq.([9])). 

The following effects will, however, lead to an overestimation of the thermal conductivity: 
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TABLE II: Thermal conductivity at T=600 K, from GK-MD compared to lattice dynamics in the 
classical limit {n{u)) — )• ksT/huj) with an equivalent number of k-mesh, and experiment for two 
different supercell sizes. 



Supercell size 


MD-GK 


LD 


experiment 


7x7x7 
10x10x10 


37 ± 10 
43 ± 12 


32.67 

47.2 


64 ± 3 

64 ± 3 



in the classical MD simulations, the number of modes is the high-temperature limit of the 
Bose-Einstein distribution, ksT/hukx which is larger than the quantum distribution. This 
leads to a heat capacity per mode of ks and therefore an overestimate of the true heat 
capacity (see also Fig. ([H])). In a finite size cell, the allowed frequencies are quantized 
and energy conservation after a 3-phonon process can never be exactly satisfied, this will 
lead to an effectively longer lifetime for phonons, and thus also overestimate n. It is not 
easy to quantify these errors except for those due to the phonon occupation numbers. It 
is therefore possible that there is a cancellation. In our case, since only two supercell sizes 
were considered, we can not do a systematic size scaling study, but overall, due to these 
cancellations the MD-GK results seem to be weakly dependent on size, in agreement with 
previous MD simulations (see for example Table I in reference [111]]). 

Here, we must point out some discrepancy between published results on Si using the SW 
potential. Using the MD-GK method, Philpot et al. and Volz et al. |8|, |9(| find a thermal 
conductivity in reasonable agreement (to within 30%) with experiments. Broido et a/. [12], on 
the other hand, have shown by solving Boltzmann equation beyond the RTA, that nsw ^ ~ 
4Kexperiment- Recently Sellan et a/.[ll[ investigated size effects in GK-MD simulations, direct 
method, and also used lattice dynamics to compute the thermal conductivity of Si from the 
SW potential. They found that kld{T = 500K) = 132W/mK, which is only 70% larger 
than the experimental value of 80 W/mK, in contrast to Broido et a/'s 12| prediction. Their 
direct method followed by scaling predicts 93 ± 18 W/mK, and their unsealed GK value for 
a 8x8x8 supercell is (231 ± 57 W/mK). 

All these results point to the subtleties involved in extracting a reliable value for the 
thermal conductivity of bulk materials, no matter what method is used. 

To investigate this discrepancy, we used our approach to extract cubic force constants 
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Band Structure and DOS for Si 
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FIG. 3: (Color online) Phonon band structure of Si using force constants up to fifth neighbor shell. 



23l |. Left: DOS , and bottom: Rescaled 



Plus signs represent experimental data of Nelin and Nilsson 
Gruneisen parameters (100 x (7 — 8)) 

from the SW potential and used LD theory to compute the corresponding thermal conduc- 
tivity. Using the same k-point mesh, in order to avoid systematic errors, in comparison to 
FP-derived force constants, we found that at 150K the thermal conductivity derived from 
SW is 80% larger than the one derived from FP-DFT calculations. 

C. Phonons, DOS and Gruneisen parameters 

In extracting the force constants, we have limited the range of the harmonic FCs to 5 
neighbor shells, and that of the cubic and quartic terms to one neighbor shell, so that MD 
simulations can be done within a reasonable time. Using the harmonic FCs, we can obtain 
the phonon spectrum. As can be seen in Fig. (|3]) the speeds of sound and most of the features 
are reproduced with very good accuracy. It is well-known that in order to reproduce the 
fiat feature in the TA modes near the X point, one must go well beyond the fifth neighbor. 
For the band structure and the density of states (DOS), the overall agreement is good, 
except for the Gruneisen parameters of the TA branch, where our calculations, which only 
include cubic force constants up to the first neighbor shell, overestimate 7(X,TA). Based on 
Klemens' formula (Eq. ([9])), one might anticipate that our model will slightly underestimate 
the lifetime of TA modes and thus their contribution in the thermal conductivity. 
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Joint Density Of States for Si 
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FIG. 4: (Color online) Top in blue is the DOS associated with two-phonon creation or annihilation 
{DOS2), and bottom in green is the DOS associated with one phonon emission or absorption 
(DOS'2 )• In red, the contribution of the matrix elements defined in Eq. (j20p are displayed. The 
peak at 500 cm^^ is the main reason for smaller lifetimes of optical modes. 

D. Phonon lifetimes and mean- free paths 

To get an idea about the relative contributions of the matrix elements, representing the 
strength of the 3-phonon interactions, versus the phase space available for these transitions, 
characterized by the two-phonon DOS, we show in Fig. (jl]) the plots of these quantities. We 
define the contribution of the matrix elements as: 

F{u) = J2 ^i^-^kx) E \V{k\,l,2)\' (20) 

kX 1,2 

From Fig. (jl]) we can note that optical phonons have a much larger weight coming from 

the matrix element \V{kX, 1, 2) p. This explains why they have such a larger relaxation rate 

compared to acoustic modes for which the matrix elements contribution is very small. The 

two-phonon DOS is representative of the phase space available for the transitions, and is 

defined as: 

DOS^iio) = J2S{uj-iOi± W2) (21) 

1,2 

From Fig. (jl]) it can be inferred that one phonon absorption or emission {DOS2) dominates 
for low frequency phonons (acoustic), while two-phonon absorption or emission {DOS2) 
dominates at high frequencies (LA and optical). 

Next, we show in Fig. (j5]) the calculated lifetimes of the 3 acoustic and optical modes 
versus frequency for a regular mesh of kpoints in the first Brillouin zone, at T=70 and 
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277K. The results depend slightly on the number of k-mesh points used for the integration 
within the FBZ. Here, we are showing results obtained with 18x18x18 mesh, which is close 
to convergence. The normal and umklapp components of the lifetimes are separated as 
1/r = 1/r^ + 1/t^ . We can note that although the lifetimes associated with normal 
processes are in l/w^, those of umklapp processes seem to scale at low frequencies like 1/u^ 
so that the former dominates at low frequencies. This is in contrast to the first-principles 
results provided by Ward and Broido 3l| where they report that the umklapp rate is in 
u'^. Even though not explicitly mentioned in their paper 32j, fits to their data with u^ was 
almost as good as the fit with u^. In the appendix, we provide a proof why in the case of Si 
the umklapp rate would behave as oj^. 
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FIG. 5: (Color online) Lifetimes of the 6 branches in Si at 277 K versus frequency on a logarithmic 
scale. Top is for normal and bottom is for umklapp processes. The quadratic dependence of the 

acoustic modes can be noticed for normal processes, while umklapp processes seems to scale as 

l/a;3. 
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From Fig. ([5]), we can notice that at low frequencies (typically below 3 THz or 100 
cm~^ where dispersions are linear), normal rates dominate while at higher frequencies and 
typically for optical modes, umklapp processes dominate transport. 

E. Thermal conductivity from lattice dynamics 

To see what is the contribution of each MFP to the total thermal conductivity, following 



the approach of Dames and Chen[33l], we have decomposed the thermal conductivity based 
on each mode and sorted each component according to their mean free paths. One can then 
define a differential thermal conductivity and the accumulated one, which is its integral: 

dn^Akx) = -Vkx^kxCvkx 



'^(A) = — E M^kx) (22) 

kX 



Nk 



The above can be plotted versus the MFP, A, seen as an independent variable. Fig. ([6]) 
shows such contribution at 277 K. Considering the extrapolated value to be 166 W/mK, 
one can notice that MFPs extend well beyond 10 microns even at room temperature. MFPs 
longer than 1 micron contribute almost to half of the total thermal conductivity! One should 
also note that the range of MFPs in Si at least, span over 5 orders of magnitude from a 
nanometer to 100 microns at room temperature. This would be larger as we go to lower 
temperatures. 

To get an acurate estimate of the thermal conductivity, one needs to extrapolate the data 
obtained from a finite number of k-mesh points, according to Eq. ( fTOl) . The extrapolated 
thermal conductivity versus temperature is p lotted in Fig. ([71) and compared to the experi- 



mental results of Glassbrenner and Slack 



29(1 . We can notice that at 



and Inyushkin et al. 

low temperatures, boundary scattering limits the experimental thermal conductivity. The 
agreement is very good in the temperature range of 100 to 500K, after which experimental 
results decay faster due to higher order phonon scatterings which are like 1/T^ or higher. 
Our resutls are within the relaxation time approximation, but one could also go beyond and 
iteratively solve Boltzmann equation as Broido et al. have done|l6l|. They have shown that 
for Si and Ge, there would be about a further 10% increase in k. 

To assess the effect of the classical approximation, which is made in classical MD simu- 
lations, we have also compared in Fig. (|8]) for a given k-point density, the classical and the 
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FIG. 6: (Color online) Cumulative contributions of phonons to the thermal conductivity at 277 
K from the 18x18x18 k-mesh data. Left is according to the wavelengths, and right is according 
to the MFPs. Both differential and cumulative thermal conductivities are shown in blue and red 
respectively. For comparison, the extrapolated (to infinite k-mesh) and experimental k are also 
shown with horizontal lines at 166 and 174 W/mK respectively. 

Thermal conductivity of bulk Si after size scaling 
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FIG. 7: (Color online) Thermal conductivity of pure Si crystal from Eq. ()17p versus temperature. 
Inset shows the extrapolation to infinitely dense kpoints. 

quantum thermal conductivities within the RTA. They are displayed with symbols on the 
lines. The quantum one is given by Eq.( 1T7|) . and the classical one uses the same expression 
in which the Bose-Einstein distribution is substituted by ksT/hu both in the heat capacity 
and in the relaxation time. We can notice that the difference is small above the Debye 
temperature, as expected, but the classical value overestimates the quantum one by 10 to 
20% as the temperature is lowered further. This is a combination of the larger heat capacity 
and a smaller lifetime in the classical case. We have also plotted the contribution of each 
mode to the thermal conductivity. We can note that at low temperatures maily the two TA 
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FIG. 8: (Color online) Quantum and classical thermal conductivities of Si versus temperature 
(from Eq. (J17p with k-mesh=18). The contribution of each mode is also being plotted. 

modes equally contribute to k, whereas at temperatures above 200 K, LA and TA modes 
equally contribute about almost 1/3 of the thermal conductivity, while LO's contribution is 
about 5%. 

The computation of the thermal conductivity using the RTA is to some extent more 
straightforward than the use of GK-MD. The former involves a double summation in the 
FBZ and has very little systematic error in it, whereas the MD simulations require an 
ensemble averaging process with a relatively large error bar, not to mention the much longer 
CPU time needed to run the MD simulations. 

For a mesh of kpoints equal to the number of primitive cells included in the MD supercell, 
we have obtained agreement between MD results and the classical version of Eq.f lT7|) . as also 
shown in table (lITll. 



VI. CONCLUSIONS 

Using first-principles calculations, we developed a classical force field which was used both 
in a molecular dynamics simulation and in the calculation of anharmonic phonon lifetimes. 
Both methods provided an estimate for the thermal conductivity of pure crystalline silicon. 
The results of these two methods agreed for the same system size in the case where k^o was 
evaluated in the classical limit. GK-MD is however much more time-consuming and includes 
large statistical errors. Furthermore it does not provide much information besides the way 
the integrated autocorrelation converges with simulation time. Size effects were discussed 
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and arguments were provided why equilibrium MD simulations converged relatively fast with 
respect to the supercell size. Lattice dynamics, on the other hand, proved to be faster, more 
accurate, and contain more useful information. The use of a linear extrapolation versus the 
inverse of the size led to a surprizingly good agreement with experiments. Such extrapolation 
is justified for relaxation rates which are quadratic in frequency at low frequencies. The 
decomposition of k, into the contribution of different mean free paths showed that in Si 
MFPs span over 5 orders of magnitude from 1 nm to 100 microns at room temperature, 
where about half of the thermal conductivity comes from MFPs larger than 1 micron. 

The developed potential has the advantage of being amenable to systematic improvement 
by including more neighbor shells at the cost of heavier calculations. The approach of using 
the FGR for the estimation of relaxation rates and the RTA or an improved approximation 
to K by solving the linearized Boltzmann equation, allows one to obtain a relatively accurate 
estimate of the thermal conductivity of an arbitrary bulk crystalline structure from a few 
force-displacement relations obtained using first-principles calculations, without any fitting 
parameters. This method paves the way for an accurate prediction of thermal properties of 
nanostructured or composite materials in a multiscale approach, which takes as input the 
relaxation times due to anharmonicity and defect scatterings. 
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VIII. APPENDIX 

In this appendix, we show the frequency-dependence of the umklapp rates. According 
to Eq. [131 the relaxation rate is a product of the 3-phonon matrix element |V^(q'A, 1, 2)p, a 
combination of occupation factors, and delta functions reflecting the constraints of energy 
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conservation. We will separately discuss the frequency dependence of the matrix element 
and the phase space term. 

First, the sum over the second momentum 2 is cancelled by the constraint of momentum 
conservation, so that the relaxation rate is just the 3D integral over qi in the FBZ. One of 
the dimensions can be integrated over by using the identity 



d?qi 6{qi - qo)/\vq,x^ - Vq+g^Xilfiqi^i) = 

J d'^SqJ/\Vq^X^ - Vq+q^X2\fiqoXl) (23) 

where qo is the solution oi Uqx + ^qoXi ~^-q-qo\2 = 0- Note that the denominator containing 
the group velocities is not small as Ai as long as A2 refer to two different branches; but in 
case Ai = A2, the denominator becomes linear in q. 

Second, for umklapp processes, in the small u limit, we must have both gi and g2 = —Q—Qi 
near the Brillouin zone boundary such that gi is inside the zone and g2 outside; so that the 
corresponding frequencies are not infinitesimally small, but their difference would be. In 
general, this forces the gi surface integral to be limited to a pocket of dimensions q located 
at the FBZ boundary, so that the surface integral is of the order of q^. But in case where 
there is a degenerate band at the zone boundary, the surface would be of order q instead. 
Different cases based on the symmetry of the crystal and the type of degeneracy have been 



discussed in detail by Herring 3J]. In our case of interest, namely Si, it is possible to have a 
3-phonon process involving a small momentum q acoustic mode connecting the LA branch to 
the LO one, with which it is degenerate, near the Brillouin zone boundary all along X ^ W, 
with a surface area Sg^, therefore, of order q. 

Third, among the two types of terms: phonon u decaying to Ui + 002 and one phonon 
absorption u + ui = UJ2, the former cannot contribute because a; ~ and ui and a;2 are 
finite. Therefore only the terms (^2 — ^i) x [6{ui — U2 + (jj) — S{ui — 002 — oj)] contribute 
to the umklapp lifetimes at small frequencies. In the latter, one can substitute rii — 77-2 by 
± u dn/duji ^ 0{q). We must remember to substitute the argument u in the relaxation rate 
by its on-shell value Wg = f x g — )■ 0. So that, in the limit of low frequencies, the inverse 
lifetime can be written as: 

d'S{q,) ujqx^\V{q, g„ -g - g„)p (24) 

\Vqo\^ - Vq+q^X2\ OUJo 
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Finally, due to the odd parity of the cubic force constants, one can show that for small q 
we have \V{q, Qo, —q — qo)\ oc Sin qRj ^fUJ^ oc y^g. 

Putting everything together, we find that the umklapp rates at low frequencies are, to 
leading order, of the form: 

^oc'occ.' (25) 

This is in agreement with our numerical findings. 

For normal processes, there is no restriction for modes 1 and 2 to be near the BZ boundary. 
For instance, in the (LA — )■ LA + TA) process the term (1 + ni + n-i) contributes and will 
not be linear in q. In such cases the rate would be in (f and would dominate terms with 
higher powers of q. 
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